scripts/Empirical example NatMath/NatMathGPCM.R

In joakimwallmark/PolyOptimalIRT:

library(mirt)

# Fit mirt model and plot -----------
titlestr  <- "National Math Test"
U         <- scan("data/NatMath.txt","o")
N         <- length(U) # Number of examinees
Umat      <- as.integer(unlist(stringr::str_split(U,"")))
n         <- length(Umat)/N # Number of items
U         <- matrix(Umat,N,n,byrow=TRUE)
U         <- data.frame(U)


natMathGpcm = mirt(data = U, model=1, itemtype = 'gpcm', SE=T)
load("data/NatMath_fittedmodel.RData")
load("data/NatMath_fittedmodel_mirt.RData")
# save(U, natMathGpcm, file="data/NatMath_fittedmodel_mirt.RData")

itemfit(natMathGpcm)

for(i in 1:length(U)){
    ItemPlot <- itemfit(natMathGpcm,
                        group.bins=15,
                        empirical.plot = i,
                        empirical.CI = .95,
                        method = 'ML')
    print(ItemPlot)
}


itemplot(natMathGpcm, item = 20, theta_lim = c(-4, 4))
temp <- fscores(natMathGpcm,
                method = "EAP",
                response.pattern = U)[, "F1"]
expect <- expected.test(natMathGpcm, matrix(temp, ncol = 1))
rowSums(U)-expect

AnalyzeResult <- Analyze(theta, thetaQnt, NatMath_dataList, ncycle=ncycle, itdisp=FALSE)
head(sort(NatMath_dataList$percntrnk))
initThetas <- percentileRankFromSumScore(U)
initThetas <- rank(rowSums(U))*100/length(rowSums(U))
head(sort(initThetas))
estimates <- thetafun(initThetas, WfdList = AnalyzeResult$parList[[10]]$WfdList,
                      dataList = NatMath_dataList)$theta_out


# do our items fit the data? High p -> worse fit
# 27 is one of the worse fitting poly, 20 is one of the best
itemfit(natMathGpcm)
itemfits2 <- itemfit(natMathGpcm,fit_stats = "X2", method = "EAP")
head(itemfits2[order(itemfits2$p.X2, decreasing = F), ])
head(itemfits2[order(itemfits2$p.X2, decreasing = T), ])
itemfits <- itemfit(natMathGpcm, method = "ML")
head(itemfits[order(itemfits$p.S_X2, decreasing = F), ])
head(itemfits[order(itemfits$p.S_X2, decreasing = T), ])

# look at them together with bins
itemfit(natMathGpcm, group.bins=22, empirical.plot = 27, method = 'ML')
itemfit(natMathGpcm, group.bins=22, empirical.plot = 20, method = 'ML')

# get center of each bin used for curve estimation
item <- 27
bincent1 <- AnalyzeResult$parList[[1]]$binctr
bincent2 <- AnalyzeResult$parList[[10]]$binctr
sbinval1 <- AnalyzeResult$parList[[1]]$WfdList[[item]]$Wbin
sbinval2 <- AnalyzeResult$parList[[10]]$WfdList[[item]]$Wbin
pbinval1 <- AnalyzeResult$parList[[1]]$WfdList[[item]]$Pbin

# bins 0:100
plotICCcombo(item, natMathGpcm, AnalyzeResult$parList[[10]]$WfdList,
             NatMath_dataList$scrfine/60*100,
             bincenters = bincent1, binvalues = pbinval1)
plotICCcombo(item, natMathGpcm, AnalyzeResult$parList[[10]]$WfdList,
             NatMath_dataList$scrfine/60*100,
             bincenters = bincent2, binvalues = sbinval2, surprisal = T)
# inf
plotICCcombo(item, natMathGpcm, AnalyzeResult$parList[[10]]$WfdList,
             NatMath_dataList$scrfine,
             bincenters = bincent1/100*60, binvalues = pbinval1, infscale = T)
plotICCcombo(item, natMathGpcm, AnalyzeResult$parList[[10]]$WfdList,
             NatMath_dataList$scrfine,
             bincenters = bincent1/100*60, binvalues = sbinval1,
             infscale = T, surprisal = T)
# 0:60
plotICCcombo(item, natMathGpcm, AnalyzeResult$parList[[10]]$WfdList,
             NatMath_dataList$scrfine,
             bincenters = bincent1/100*60, binvalues = pbinval1)
plotICCcombo(item, natMathGpcm, AnalyzeResult$parList[[10]]$WfdList,
             NatMath_dataList$scrfine,
             bincenters = bincent1/100*60, binvalues = sbinval1, surprisal = T)